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ABSTRACT 


The objective of this fluids analysis is to model propellant slosh for the Europa Clipper mission 
using a two-pendulum model, such that controls engineers can predict slosh behavior during 
the mission. Propellant slosh causes shifts in center of mass and exerts forces and torques on 
the spacecraft which, if not adequately controlled, can lead to mission failure. The two- 
pendulum model provides a computationally simple model that can be used to predict slosh for 
the Europa Clipper tank geometry. The Europa Clipper tank is cylindrical with a domed top and 
bottom and includes a propellant management device (PMD). Due to the lack of experimental 
data in low gravity environments, computational fluid dynamics (CFD) simulation results were 
used as “real” slosh behavior for two propellants at three fill fractions. Key pendulum 
parameters were derived that allow the pendulum model’s center of mass, forces, and 
moments to closely match the CFD data. The parameter trends were examined as a function of 
tank fill fraction and compared with solutions to analytic equations that describe the frequency 
of slosh in tanks with simple geometries. The trends were monotonic as expected, and 
parameters resembled analytical predictions; any differences could be explained by the specific 
differences in the geometry of the tank. This paper summarizes the new method developed at 
Goddard Space Flight Center (GSFC) for deriving pendulum parameters for two-pendulum 
equivalent sloshing models. It presents the results of this method and discusses the validity of 
the results. This analysis is at a completed stage and will be applied in the immediate future to 
the evolving tank geometry as Europa Clipper moves past its preliminary design review (PDR) 
phase. 


INTRODUCTION 


When liquid propellant sloshes in a tank, it shifts the center of mass (CM) of spacecraft and 
exerts forces and torques on the spacecraft. This motion must be predicted and accounted for 
in the spacecraft design to achieve mission success. In the past, NASA’s Solar Dynamics 
Observatory (SDO) was forced to enter safe mode when the attitude control system was unable 
to accurately point the spacecraft due to fuel slosh (Mason, 2011). If slosh motion can be 
accurately predicted, attitude control systems (ACS) thrusters can be sized to adequately 
counteract the forces and torques and resonant frequencies of the system can be avoided. The 
ability to predict slosh behavior is especially important in Europa Clipper because the propellant 
makes up such a large percentage of the total spacecraft mass. The most accurate means of 
determining slosh behavior, short of performing experiments in low gravity, is using 
computational fluid dynamics (CFD). However, these simulations are very computationally 
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expensive and can take several days to generate a few minutes’ worth of data. As a result, 
pendulum models are used as an approximation of the slosh. A few representative CFD cases 
are run to calibrate key parameters in the pendulum model, and interpolation can be used to 
approximate the pendulum parameters between CFD cases. Single-pendulum models are 
commonly used in slosh modeling, but they are less accurate than two-pendulum models for 
tanks containing a sponge propellant management device (PMD), like the ones on Europa 
Clipper. This paper summarizes a new method for determining pendulum parameters for a two- 
pendulum model. It then presents the results, examines the trends in comparison with existing 
literature, and discusses of validity of the pendulum models for one iteration of Europa 
Clipper’s propellant tanks. 


THEORY 


Computational Fluid Dynamics (CFD) Model Setup 

The CFD program STAR-CCM+ was used for the Europa slosh analysis. The CFD provided the 
center of mass, forces, and moment on the tank during the sloshing event. In addition, the CFD 
program outputted the moment of inertia of the settled propellant. This program has been 
used successfully in the past by NASA Goddard Space Flight Center to determine slosh for the 
ICESat-2 and OSIRIS-Rex missions. 


CFD simulations were run for each of the two propellants, nitrogen tetroxide (NTO) and 
monomethyl hydrazine (MMH) at three different tank fill fractions. The propellant properties 
were calculated from equations provided in the United States Air Force handbooks on 
hydrazine fuels (Marsh, 1970) and nitrogen tetroxide oxidizers (Wright, 1977). 


The liquid surface was initialized at a 15-degree angle from horizontal for all cases unless 
otherwise noted. This pendulum initial offset angle was assumed to be equal to the fluid 
surface angle. Figure 1 below shows what an offset of the liquid surface looked like. The PMD 
vanes are also shown outlined in black. 
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Solution Time O (s) 


V Volume Fraction of NTO 
zx 0.00000 0.20000 0.40000 0.60000 0.80000 1.0000 


Figure 1. Offset of liquid of 15 deg with no surface tension model turned on at a fill fraction of 
35%. 


The acceleration was applied in the CFD model using the gravity model. 


The volume of the tank was discretized using a polyhedral mesh with growth cells at the solid 
boundaries (called prism mesh cells in STAR-CCM+). A mesh independence study was done to 
determine appropriate mesh refinement and parameter selection. 


The CFD program output tables for the propellant center of mass location, forces, and torques 
exerted on the tank in three cartesian coordinates — one along the axis of the tank and two 
transverse orthogonal directions. The program also outputted the moment of inertia of the 
liquid and gas as well as images of the volume fractions of liquid and gas. 


The models were released from the initial offset and allowed to settle for sufficient time to 
determine the sloshing frequency and damping. 


Pendulum Model Derivation 


A damped pendulum’s movement can be described by the equations below. 
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Equation 1 


2 


da 
mi? ae 2fwml2 ae mglsin@ = 0 


Equation 2 
Equation 1 shows that the moment of inertia of the pendulum, J, multiplied by the angular 
_ a6. : 
acceleration, q's equal to the torque, Tt. For a damped pendulum, Equation 1 becomes 


Equation 2 where m is the mass of the pendulum, l is the length of the pendulum, g is the 
spacecraft acceleration, w is the undamped frequency of the pendulum, ¢ is the damping ratio, 
and @ is the angle of the pendulum from the resting position. The forces used to obtain 
Equation 2 are shown in Figure 2. 
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Figure 2. Forces on pendulum bob. 


A positive rotation is counter-clockwise. By dividing by ml?, by noting that w? is equal to 7 and 


by assuming a small angle approximation, Equation 3 can be obtained. 
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Equation 3 


Equation 3 is an ordinary differential equation and there are known methods for solving it. The 
result is shown below in Equation 4, with 09 being the initial pendulum offset angle. 


A(t) = Ooe-Set(— 2 
@W 


perc 


sin (w 1% t) +cos(w 1-@t)) 


Equation 4 


Note that Equation 4 predicts no slosh if 8, is equal to zero. This is why the CFD simulation 
requires the surface to be offset to one side. 


The unknowns in Equation 4 are ¢, w, and 0p. If w is determined, the length of the pendulum, l, 
can be calculated using the relation in Equation 5. 


Equation 5 


In addition, the hinge location, Ay, the static mass, mp, and the static mass location, hy need to 
be found to accurately calculate the force in the axial direction and the moments applied to the 
tank. All quantities are fixed within the model for a given propellant fill fraction, except for the 
pendulum angle, which is a function of time. 


To solve for the unknowns, the pendulum parameters are solved so the pendulum model 
output fits the CFD model outputs. 


The forces that the pendulum exerts on the tank walls are calculated in Equation 6 and 
Equation 7. Equation 6 is the summation of the forces in the transverse direction due to the 
pendulum and Equation 7 is used to calculate the summation of the forces in the axial direction 
due to the pendulum. (Note: Equation 7 can be modified to include the term —mgq, the inertial 
force on the static mass, if matching with the forces in the axial direction from the CFD. But this 
added force will not cause a moment because the static mass stays on the axial line of the tank, 
and so should not be used in Equation 8). The summation of the moments was found using 
Equation 8, where g is the spacecraft acceleration, 1, and 1, are the moment arms in the axial 
and lateral directions, respectively. In this coordinate system, the z coordinate is in the 
direction orthogonal to both a and t. 
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dt? 
Equation 6 
F, = —ml ey cos(@) — nie aes —mg 
dt dt? 
Equation 7 
M, = —F,l,+ hl 
Equation 8 


Figure 3 illustrates the forces exerted by the pendulum on the tank. 
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Figure 3. Forces exerted by pendulum on tank. 


The parameters m, ¢, and w (components of @(t)) were selected such that the pendulum 
motion closely matches the CFD data. 


The hinge point, h,, for each pendulum in the model is chosen so that M, from the pendulum 
model approximately matches the M, from the CFD model. The CFD model sums the moments 


about z at a=0 and t=0. The CFD model’s origin is located at the top of the tank. The pendulum 
model uses the same coordinate system. 
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The magnitude of the static mass, mp, and the static mass location, ho, are calculated using 
Equation 9 and Equation 10, respectively. The center of mass location of the fully-settled 
propellant, h;, and total mass, m7, are taken from the CFD model results. 


Mj = ™Mr-m 
Equation 9 


mrhr — mh 
No TT p 


Mo 


Equation 10 


The equations can be modified to find results for a multiple pendulum model. To find the 
forces and moments, the forces and moments from each individual pendulum are added 
together. To find the static mass magnitude, m is replaced with (m, + m,...+m,) in 
Equation 9. To find the static mass location, mh, is replaced with (m,hpy; + Mzhpp ... + 
MnNpn) in Equation 10. 


The static mass is treated as both a point mass and a body with a moment of inertia about its 
center of mass. This approach is necessary so that the pendulum model’s total moment of 
inertia is equal to the propellant’s moment of inertia when the propellant is in a settled 1-g 
configuration, [n7-op. While during sloshing the moment of inertia of the body will change due 
to the change of shape due to sloshing, the settled 1-g propellant configuration is an adequate 
approximation. 


= 2 
I prop ae IpropcM ote MrD5rop 


Equation 11 


Ipropcm is the moment of inertia of the propellant about its center of mass, outputted by CFD. 
Dyrop is the distance between the propellant center of mass and the static mass location. 
Equation 11 is the parallel axis theorem, used to calculate the moment of inertia of the 
propellant mass about the static mass location. The pendulum mass moments of inertia, [j,, 
are calculated about the static mass location using Equation 12, where D,,, is the distance 
between the pendulum mass and the static mass. 


bo=-mDe2, 


Equation 12 


Using a two-pendulum model, Equation 13 shows the equation for calculating Imoem, the 
moment of inertia of the static mass about its own center of mass, where the subscripts 1 and 2 
refer to the two pendulums, respectively. 
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= 2 2 
Imocm — Iprop ~ ™mDmi1 _ M2zDmo2 


Equation 13 


Since the pendulums, when resting, are aligned with the center-line of the tank, the moment of 
inertia of the pendulums are zero when calculated about the center-line of the tank. 


Basis for Two-Pendulum Model 


NASA’s Cassini mission used PMD propellant tanks similar to those on Europa Clipper, and also 
used a two-pendulum model for slosh analysis. A notional tank and PMD is shown in Figure 4. 


Figure 4. Notional tank and PMD 


In a paper written by Enright and Wong on Cassini slosh analysis (Enright, 1994), the authors 
used two separate pendulums to model two modes of propellant sloshing, which they called 
“sector” and “full” tank slosh modes. They base this on the observation that the tank geometry 
produced by the PMD can be approximated as both a sector of a cylinder and an annular 
cylinder, as shown in Figure 5. The two pendulums in the model account for the two slosh 
modes between the PMD vanes and around the PMD vanes. 


Sector tank mode (top view) Annular tank mode (top view) 


PMD 


Tank Wall —¥ 
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Figure 5. Schematic of PMD approximated as annular and sector tank modes. 


Extensive analytical and experimental work was conducted on liquid sloshing in simple 
geometry tanks including cylindrical sector and annular tanks. Much of this work was compiled 
in NASA’s technical report SP-106 (Abramson, 1966), and updates were later added by the 
Southwest Research Institute (SwRI) (Dodge, 2000). These documents have been used 
extensively in the past by slosh analysts. They documents provide analytical equations for 
frequency associated with each of these slosh modes, which are used later as a sanity check to 
validate the model parameters. Additionally, equations from these documents for damping 
ratio in a bare cylindrical tank were used as a comparison for model parameters. 


METHODOLOGY 


High Level Summary 

The goal was to find the parameters of two damped pendulums and one static mass that, when 
superimposed to calculate a total CM, most closely matched the CFD CM over time. This was 
done for two propellants at three fill fractions each, with the goal of selecting parameters that 
were monotonic as a function of fill fraction to allow for interpolation between fill fractions. 
The behavior before the first peak (the first wave of movement from one side of the tank to the 
other) experienced highly nonlinear damping, so a pendulum model was not able to capture 
the behavior before the first peak while also accurately modeling the behavior after this first 
peak. Therefore, the decision was made to fit the pendulum model to the behavior after the 
first peak only. It will be shown later that this is a conservative choice for all initial offset angles 
greater than the 15-degree offset used in the CFD simulations. 


Examination of the CFD CM data reveals that it can be approximately resolved into the sum of 
two separate damped oscillations. Initial guesses for each pendulum’s frequency w and 
damping ratio ¢ were calculated by selecting points in the CFD data that were approximately 
the peak-to-peak values for each pendulum. A mass fraction mf was defined to describe the 
pendulum bob mass as a fraction of the total propellant mass. Then, the first pendulum’s mfz, 
W@1, and (1 were generated using curve fitting to the CFD CM data. The second pendulum’s mfz, 
W2, and ¢2 were generated by curve fitting to the difference between CFD CM data and the first 
pendulum’s CM contribution. Additional steps were taken to fine-tune the parameters to 
improve the fit and minimize the error in force. Pendulum attachment height h was found by 
minimizing the error in moment. Finally, moments of inertia /, and /, were assigned to the static 
mass to have the model reflect the total moment of inertia of the settled fluid. 


Initial Guesses 


Figure 6 shows an example of the CFD CM data and the peaks used to find initial guesses for w 
(used along with acceleration to calculate /ength) and Z. 
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Center of mass in x vs. 


T T T T 


Figure 6. CFD CM data with peaks indicated. 


The first set of consecutive peaks is selected at the beginning of the simulation where the 
longer pendulum’s movement dominates the CM movement. Equation 14 and Equation 15 are 
used to calculate w, and @: from these points. The second set of consecutive peaks is selected 
near the end of the simulation, after the longer pendulum has damped out and the shorter 
pendulum’s movement dominates the CM movement. wz and @ are calculated from these 
points. In some data sets, the CFD CM oscillates about a non-zero center. This offset is on the 
order of 1 mm or less and is likely due to numerical error in the simulation. To center these data 
sets about zero for the initial guess calculation, the average CM of the last few seconds was 
subtracted from every point in the CM data. 


Equation 14 


2 
a ees 
In(CM,/CM2) 
Equation 15 


Ti and 72 are the times of the first and second consecutive peaks, respectively, CW1 and CM2 
are the center of masses of the first and second consecutive peaks, respectively, and 6 is the 
logarithmic decrement. 
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Curve Fitting Pendulum Equation 

Matlab’s curve-fitting fsolve function was used to iteratively solve for a more accurate mass 
fraction mf, frequency w, and damping ratio ¢ that would allow the pendulum model to fit the 
CFD data. The initial guesses passed into fsolve were 0.5 for mf, and the previously generated 
initial guesses for w and ¢. The equation passed to fsolve was Equation 16, where m was the 
mass of the pendulum bob, m: was the total mass of the propellant, L was the length of the 
pendulum, and @ was the angle of the pendulum over time given in Equation 4. 


m mm. 
CMx(t) = 7 EM) = me sin(@(t)) 
t t 
Equation 16 


fsolve was used first to fit a longer pendulum to the high-amplitude oscillation, as shown in 
Figure 7. The blue circles indicate the points of interest where fsolve aimed to match the 
pendulum’s CM contribution with the CFD CM. The points were selected to span the time 
starting at the first peak and ending when the oscillation appeared to damp out completely. 


Center of mass in x coordinate 
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Figure 7. Long pendulum fit to CFD CM data. 


Then, fsolve was used to match a second pendulum to the difference between the first 
pendulum’s CM contribution and the CFD CM. The difference is shown in black in Figure 8. The 
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points of interest here were selected in a time interval where the difference is close to the 
shape of a damped sinusoid. 


Center of mass in x coordinate 
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Figure 8. Short pendulum fit to the difference between CFD CM and the long pendulum’s CM 
component. 


For both pendulums, trial and error was required to find the points of interest that allowed 
fsolve to produce an accurate fit. This included changing the time interval or the timestep 
between points. For the shorter pendulum especially, it was often necessary to subtract a 
correction factor to center the CFD oscillations about zero. 


After mf, w, and ¢ were determined for both pendulums, a total CM was calculated by taking 
the sum of Equation 16 for both pendulums. 


Iterations to find mass participation fractions 

For some test cases, it was possible to further improve the pendulum fit by iterating on the 
mass fraction parameter. The underlying reasoning is that if the model were perfect, then when 
one pendulum was passing through vertical, the other pendulum’s CM contribution would 
account for the entire CFD CM of the fluid. Using this reasoning, several points of interest were 
selected where the first pendulum was vertical. This happens when the first pendulum’s CM 
contribution crosses zero (intercepting the time axis). In Figure 9, the black x’s are points where 
Pendulum 2’s CM crosses zero; so if the model were perfect, Pendulum 1 would equal CFD. 
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CM in X vs. Time 
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Figure 9. Points of interest for iterating on mass fraction. 


Then, Equation 17 was used at each of these points to calculate the mass fraction required of 
the second pendulum to match the CFD CM. These mass fractions were averaged and used as 
the new parameter for the next iteration. The same calculation was carried out for points of 
interest where the CM of the second pendulum crossed zero. Iterations were performed until 
the mass fractions converged or it became clear that they would not converge. 


CMx 


Mi Lz sin 62 


Equation 17 


Additional Fine-Tuning 

After using curve fitting to obtain a full set of parameters for each propellant at each fill 
fraction, the parameters were fine-tuned to improve the fit and meet certain criteria. The 
following adjustments were only used in cases where minor adjustments could achieve the 
criteria without compromising the quality of fit. The values for m, w, and ¢ were adjusted to be 
monotonic as a function of fill fraction. The values for m and w were adjusted to be identical 
between the two different propellants. The plot of the model vs. the CFD data was examined 
closely to adjust the parameters to fit the peaks more accurately. The final model for 50% fill 
fraction is shown in Figure 10. 
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CM in X vs Time 
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Figure 10. Pendulum model vs. CFD CM at 50% fill fraction. 


Quantifying Accuracy of Fit 


The metric used to quantify the accuracy of the model was the mean difference between 
pendulum force and CFD force. This was calculated using Equation 18, where n is the number of 


points in the range of interest. The black arrows in Figure 11 show the force errors that were 
averaged. 
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n 


1 
5 ». abs(CFD — pendulum) 


1 


Equation 18 


Force in X vs Time 
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Figure 11. Illustration of error in force. 


It is important to note that a slight time shift between the CFD peaks and pendulum model 
peaks can lead to a misleadingly high error. A visual inspection of the matching between the 
CFD and pendulum outputs was done as an added check to ensure a good fit. 


Finding Pendulum Attachment Point 

The “hinge height” is the point along the tank’s central axis where the pendulum string should 
be attached to most accurately model the moment exerted by the fluid. Equation 8 is rewritten 
as Equation 19 to include hinge height h as a variable. 
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M, = Fla + Fil, 
= (ml 6? sin@ — ml 6 cos @)(—(h — Lcos @)) + (—ml6? cos 6 — ml6 sin @ 
—mg)(lsin 0) 


Equation 19 
M. 


M, Z41 


+ M,, 


total = 


Equation 20 


Initially, an attempt was made to use fsolve to fit a curve to the CFD moment data while 
allowing hi and hz2to vary. However, this method did not result in a good fit. Instead, hz was 
varied by manually inputting different values. The value was selected that minimized the mean 
error in moment (calculated using the same method as mean error in force in the previous 
section). This process was then repeated for hz. In some cases, the optimized hinge heights 
resulted in the pendulum bob being outside the physical confines of the tank walls, and the 
hinge height was adjusted slightly to avoid this, but only if quality of fit was not significantly 
impacted. An example of the pendulum bob inside the tank walls is shown in Figure 12. 
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Figure 12. Pendulums at maximum displacement within tank walls. 


Finding Static Mass Moment of Inertia 

In order to accurately model the total moment of inertia of the propellant, a moment of inertia 
is assigned to the static mass, as described in the theory section. This calculation used the 
settled state of the propellant and pendulums and the reference point was the center of mass 
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of the settled propellant. The static mass’s moment of inertia about the axial axis is equivalent 
to the settled propellant’s moment of inertia about the axial axis (calculated using CFD). The 
static mass’s moment of inertia about the z-axis is given by solving for [static in Equation 21. 
Again, the z-axis is the direction orthogonal to both the tank’s vertical axis and the lateral 
sloshing axis. 


Isettled propellant about propellant CM — In settled" pendulm model about propellant CM 
2 Z 
=™, (CH, —l)- CM prop) + M2 (CH, =) = COM otop) 
2 
on (Irieatie + ™Mo (Hstatic ~~ CM ian) ) 


Equation 21 


Figure 13 illustrates the variables in Equation 21. 


Figure 13. Schematic of tank with settled pendulums. 


RESULTS AND DISCUSSION 


Results 


Table 1 below outlines each of the CFD cases that were run. Cases 1-3 represent different fill 
fractions of NTO and Cases 2-6 represent different fill fractions of MMH. The 2_8 cases are the 
test cases in which Test 2’s CFD surface offset angle is changed to 8. These cases are not used 
in determining trends as a function of fill fraction, but give essential insight into the limits of the 
pendulum models. 
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Table 1. High-g Slosh Test Cases for November 2016 Europa Slosh Deliverable 


Case Acceleration | Temperature | Initial Fill Mass (kg) | Propellant 
(m/s) (deg C) Offset Fraction 

Angle (%) 

(deg) 
1 0.071 15 15 25 418.6 NTO 
2 0.057 0 15 50 855.5 NTO 
3 0.039 0 15 85 1454.3 NTO 
4 0.071 15 15 25 252.6 MMH 
5 0.057 0 15 50 513.3 MMH 
6 0.039 0 15 85 872.6 MMH 
2_10 0.057 0 10 50 855.5 NTO 
2_5 0.057 0 5 50 855.5 NTO 
2_30 0.057 0 30 50 855.5 NTO 


The pendulum parameters derived for these cases are listed below in Table 2. These 
parameters prioritized fitting the CFD slosh behavior after the first CM peak and ignored the 
behavior before the first peak. The Methodology section discusses the reasoning for fitting only 
one of these regions. 


Table 2. Pendulum Model Parameters for Different Test Cases, Fitting After First Peak 


Test Case Number 
1 2 3 4 5 6 
Propellant NTO MMH 
Fill Fraction (%) | 25 50 85 25 50 80 
Accel (m/s?) 0.071 0.057 0.039 0.071 0.057 0.039 
[mi(kg) | 20.09 [4449 [21087 |1212 [2669 | 126.53 | 
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Mz (kg) 12.56 24.81 26.18 7.58 14.89 15.71 
Mo (kg) 385.95 786.20 1217.25 232.90 471.72 730.3662 
mf, (kg) 0.048 0.052 0.145 0.048 0.052 0.145 
mfz (kg) 0.03 0.029 0.018 0.03 0.029 0.018 
mfp (kg) 0.922 0.919 0.837 0.922 0.919 0.837 
hy (m) 0.9 0.4 0.5 0.9 65 0.5 

hy (m) -1.0 0:7 0.3 0.9 0.7 0.2 

ho (m) -1.12 -0.99 -0.79 -1.14 -0.99 -0.8 

L, (m) 2.118 0.651 0.353 2.118 0.651 0.353 
L, (m) 0.140 0.132 0.301 0.140 0.132 0.301 
& 0.340 0.105 0.035 0.350 0.110 0.037 
& 0.015 0.022 0.035 0.020 0.025 0.037 
fi; (Hz) 0.029 0.047 0.053 0.029 0.047 0.053 
fa (Hz) 0.113 0.105 0.057 0.113 0.105 0.057 
I7-static mass 

(kg-m?) 32.52 88.89 224.06 19.82 52.88 133.75 
Iq—static mass 

(kg-m?) 55.97 128.54 228.64 33.87 77.17 137.20 


There is not a unique solution when curve fitting to the results. Parameters were chosen that 
provided a good fit, but since the pendulum model cannot match the results exactly, 


engineering judgment was used to determine what constituted a good fit. The propellant slosh 


before the first CM peak was ignored during the fitting process and thus is not reflected in the 


pendulum model. 


Figure 14 through Figure 22 illustrate the quality of fit of the pendulum model for some of the 


test cases. 
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CM in transverse direction vs. time for 25% fill fraction. 
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Figure 15. Force in transverse direction vs. time for 25% fill fraction. 
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Figure 16. Moment about z-axis vs. time for 25% fill fraction. 


—— CFD 
Pendulums 

——- pendulum 1 

— —~pendulum 2 


CM(m) 


Time(s) 


Figure 17. CM in transverse direction vs. time for 50% fill fraction. 
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Figure 18. Force in transverse direction vs. time for 50% fill fraction. 
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Figure 19. Moment about z-axis vs. time for 50% fill fraction. 
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Figure 20. CM in transverse direction vs. time for 80% fill fraction. 


5 ——cFD 
Pendulums 
4 ——-pendulum 1 
— —=pendulum 2 
3 1 
| 
| 
—~ 2h 
= ik \ 
8 | Y fA 
1 p 
f /’ A A P 
y \)\ y 
0 7 ‘oy ~ y ~— = \ ~ ~y= i f_ is \ 
\\ / \ {/ al 
=4 
VY 
-2 
0 20 40 60 80 100 120 140 160 


Time(s) 


Figure 21. Force in transverse direction vs. time for 80% fill fraction. 
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Figure 22. Moment about z-axis vs. time for 80% fill fraction. 


Discussion of Trends and Literature Comparison 


The pendulums’ mass, frequency, damping ratio, and hinge height were plotted as a function of 
fill fraction to examine the trends. These trends were compared to analytical and experimental 
trends for simpler tank geometries presented by Abramson (1966), Enright (1994), and Dodge 
(2000). 


Figure 23 shows that the pendulum bob mass fractions are monotonic as a function of fill 
fraction. It was found that identical mass fractions for NTO and MMH at a given fill fraction 
were able to accurately model both propellants. The trend indicates that for interpolating 
between fill fractions, a piecewise linear fit would be more accurate than a linear fit. The 
physical reasoning for this is that for the first two fill fractions, the PMD is partially covered by 
the propellant (see Figure 1), so sloshing occurs primarily between the vanes. For the last fill 
fraction, the PMD is completely covered, so a different slosh behavior occurs above the PMD. 
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Figure 23. Mass fraction as a function of fill fraction. 


Figure 24 plots pendulum frequencies as a function of fill fraction, which are also monotonic as 
a function of fill fraction. As with mass fraction, it was found that identical frequencies for NTO 
and MMH at a given fill fraction was able to accurately model both propellants. Also shown in 
Figure 24 are analytical frequencies in a “bare tank” (cylindrical flat-bottomed tank with no 
PMD) given by the equation in Section 3.2 of Dodge (2000). A comparison of the two pendulum 
frequencies to the bare tank frequency was also included for a slightly different tank geometry 
in the Cassini mission (Enright, 1994), shown here as Figure 25. Both plots show the two 
pendulum frequencies converging to the bare tank frequency as fill fraction increases. The 
physical explanation for this is that as the fill fraction increases, the PMD becomes fully covered 
by the propellant, so the slosh behavior begins to resemble that of a bare tank with no PMD. 


TFAWS 2017 — August 21-25, 2017 25 


Frequencies vs. Fill Fraction 
Pendulum Fits vs. Analytical Bare Tank Equation 
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Figure 24. Frequency vs. fill fraction compared to analytical bare tank. 
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Figure 25. Cassini frequency as a function of fill fraction (Enright, 1994). 


Abramson (1966) provided analytical equations for sector and annular frequencies in a 
cylindrical tank in SP-106. Figure 2.3 in Section 2.3 of SP-106 shows how h, a, b, and a are 
derived from the geometry of the sector or annular tank. Equation 2.19 in Section 2.3 of SP-106 
shows the equation for slosh frequency as a function of acceleration, fluid height, and tank 
geometry, copied here as Error! Reference source not found.. g is the acceleration. The 
coefficient (mn in this equation is tabulated by Bauer (1963) and is dependent on tank 
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geometry. The subscript m is 1 for the mode of interest, as explained by Dodge (2000) Section 


1.5. 


Equation 22 


Table 3 summarizes the variables used in Equation 22 to calculate frequencies. 


Table 3. Parameters Used in Analytical Sector and Annular Frequencies 


Parameter | Sector Slosh Annular Slosh 

a (m) Radius of tank Radius of tank 

b(m) 0 Average radius of PMD vane 

a 1/8 1/2 

P oe 2 

k b/a b/a 

C10 Look up from Bauer (1963) tables | Look up from Bauer (1963) tables 


The analytical sector and annular frequencies calculated using this equation were plotted in 
Figure 26 and compared to the pendulum model frequencies. The plot shows that the 
pendulum frequencies and analytical frequencies are close. The physical reason for the 
differences include the fact that the PMD is neither a true sector nor annular tank, which will 
affect all fill fractions. The inner radius of the annular tank was approximated as an average 
radius of the PMD vane, but the PMD shape is not physically representative of an annular tank. 
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Frequencies vs. Fill Fraction 
Pendulum Fits vs. Analytical Sector and Annular Tanks 
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Figure 26. Frequencies as a function of fill fraction compared to analytical sector and annular 
tanks. 


There was a larger discrepancy in pendulum and analytical annular slosh frequencies at 25% fill 
(circled in yellow). This was likely due to the fact that the CFD data came from a domed bottom 
tank, while the analytical equation was for a flat-bottomed tank. At 25% fill, the propellant is 
almost entirely sloshing in the domed section, so that geometry difference will affect the 
frequencies more than at higher fill fractions. There was also a larger discrepancy in the sector 
slosh at 85% fill (circled in green). This was likely due to the PMD being fully covered by the 
propellant at this fill fraction. The sloshing occurred above the PMD where the tank geometry is 
similar to a bare cylinder without vanes to form sectors. 


Damping ratio trends were also validated against correlations given in literature. Dodge (2000) 
Section 2.2 includes a correlation equation for damping ratio in a bare cylindrical tank (no 
sectors or annulus) as a function of geometry, fluid height, and acceleration. It was developed 
by Mikishev and Dorozhkin based on their experimental tests and copied here as Equation 23. It 
includes a correction coefficient Cgome to correct for a dome-bottomed tank, which can be read 
off a chart labeled Figure 2.2 in Dodge (2000), copied here as Figure 27. For the 25% fill fraction 
Cdome = 2 was used, and for 50% and 85% fill fraction, Caome = 1 was used. These values for Cgome 
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correspond to liquid fill ratios of h/R, where h is the depth of the settled propellant from the 
end of the hemispherical bottom and R is the radius of the cylindrical portion of the tank. 


h 
= Caome0-79 Re, |1 oes _ 7k 
Y = Caome- ot ae » (Ea) * n (E84) 
sinh (—3— cosh {—3— 
Equation 23 
= 8.0 
ron 
So 
=. 66 
E 
oO 
z= | 
® 4.0 
§ 
S 
= 2.0 
G 
S 
S OOO 0.4 1.6 2.0 


0.8 1.2 
Liquid fill ratio, h/R 
Figure 27. Correction coefficient for domed bottom tank. 


Figure 28 shows damping ratios for each pendulum and the damping ratios given by Mikishev 
and Dorozhkin’s correlation. The damping ratios of each pendulum were monotonic as a 
function of fill fraction. The more viscous propellant, MMH, has a slightly higher damping ratio 
than NTO. The damping ratios for two pendulums are somewhat close to the damping ratio in a 
bare tank, especially as the PMD becomes covered and the slosh behavior becomes more like 
sloshing in a bare tank. 


The pendulum model’s increase in damping due to increased viscosity (i.e. increase from NTO 
to MMH) is much less than what the analytical correlations would predict. This could be related 
to the fact that the damping caused by a PMD is dominated by a drag force, while the damping 
in a bare tank is dominated by viscous force. Since the PMDs are identical for the two 
propellants, the drag force is identical, and the change in viscosity has only a small impact on 
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the damping coefficient. This idea was based on Section 2.3 of Dodge (2000), in a discussion of 
damping by ring baffles. 


Damping Ratio vs. Fill Fraction 
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Figure 28. Damping ratio as a function of fill fraction compared to analytical. 


Dodge (2000) presents an additional correlation by Stephens, et al. which yields similar bare 
tank damping ratios, but differs in that it doesn’t provide a correction coefficient for the domed 
bottom. There was no analytical equation found for damping ratio in sector and annular tanks. 


The final parameter trend is the height of the pendulum attachment point along the axis of the 
tank, referred to as “hinge height” and plotted in Figure 29. The values are measured from a 
reference point at the top of the tank. 
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Hinge height vs. fill fraction 
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Figure 29. Hinge and static mass height as a function of fill fraction. 


The pendulum length was calculated from the frequency and acceleration using the relation 


oe 7 These values were used in conjunction with hinge heights to draw a schematic of the 


pendulums within the tank, as in Figure 30. From this schematic, the 50% and 85% fill fraction 
pendulums were all within the region physically occupied by the propellant. In the first iteration 
of hinge heights, Pendulum 1 for 25% fill was below the tank, so the hinge height was shifted up 
slightly to keep it within the tank. This shift had minimal impact on the accuracy of modeling 
moments, since the moments are highly insensitive to changes in hinge height. 
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Figure 30. Schematic of pendulums within tank walls. 


Discussion of Validity of Assumptions 


The analysis in this paper made the following assumptions: 


The initial pendulum offset angle is equal to the CFD liquid surface offset angle. 


The assumption that the initial pendulum offset angle is equal to the CFD liquid surface offset 
angle may be incorrect, as will be shown shortly, but there is no other indication of what the 
initial pendulum offset angle should be. Several test cases were run in which the initial 
pendulum offset angle was changed to verify if the specific initial pendulum offset angle would 
give a “best” result. The pendulum parameters did vary for different initial pendulum offset 
angles, but there was no clear “best” result. To better understand the physics that might allow 
the pendulum data to match with such a wide variety of pendulum parameters, several more 
CFD cases were run with different surface initial angle offsets. The results are shown below in 
Figure 31. 
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Figure 31. Comparison of CFD forces in the x-direction for different initial surface offset 
angles. 


The CFD results from the 5 deg through 15 deg offset angle cases show that after an initial spike 
in the force, that the slosh forces become almost identical; however, at an initial offset angle of 
30, the behavior is significantly different. It takes several periods before the influence of the 
initial spike is dampened out for the 30 deg case, with a phase shift as the only remaining 
difference. Several physical explanations could be put forward, though more research would 
be needed to prove what is going on. One possible explanation is that the initial sloshing event 
involves a large percentage of the propellant moving to the bottom of the tank and the majority 
of this propellant motion is quickly dampened out by the vanes in the PMD. After that, a 
smaller bulk slosh is left along with sector slosh. The damping of the remaining pendulum slosh 
is drastically different from the damping in the initial slosh event because less propellant is 
interacting with the PMD. Additionally, fluid damping is a function of wave height, and the 
initial condition is a large wave that dampens out much more quickly than the small subsequent 
waves. As the initial surface offset angle is increased, the influence of the sloshing waves 
cresting and then rejoining the bulk liquid, as well as the inability of the PMD to fully dampen 
out the initial bulk movement quickly, leads to non-linearities that cannot be represented by a 
linear pendulum model. The results from the cases at 25% and 50% fill fraction show this initial 
spike, but for 80% fill where the PMD is fully immersed, this initial spike does not seem to be 
present. A large number of CFD slosh cases would need to be run over a range of fill fractions, 
accelerations, and initial surface offset angles (as well as zero-g initial offsets) to determine how 
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the initial pendulum angle should be selected. Also, it would be necessary to make the 
damping non-linear. 


With the current data available, choosing the initial angle offset of the pendulum to equal the 
CFD liquid surface initial angle offset is the most reasonable approach since there is no other 
physical data within the CFD setup that would indicate another angle. 


In light of this, it is critical to realize that using the pendulum model with a smaller initial fluid 
angle would under-predict the CM and force amplitude in the propellant. This is because the 
pendulum parameters were derived with an initial propellant and pendulum angle of 15 
degrees, and the pendulums’ oscillation amplitude is proportional to the initial angle, while the 
propellant’s oscillation amplitude appears to be constant for initial angles between 5-15 
degrees. However, for initial disturbance angles greater than 15 degrees, the pendulum model 
should overpredict the liquid’s oscillation amplitude. 


The pendulum masses, lengths, damping, and hinge points can be determined by choosing 
them so that the pendulum results match the CFD results. 


The assumption is used because there is not sufficient theory to fix any of these values 
independent of the CFD results. Because no better theory exists to the author’s knowledge, a 
trial and error method was used. During the trial and error process it was noted that there was 
not one unique solution. If the mass was increased, the CM oscillation amplitude was higher, 
but that could be offset by increasing the damping ratio. Because of this, confidence in the 
accuracy of the mass and damping parameters is lower than the confidence in the accuracy of 
frequency. Some attempt was made to try and have the initial center of mass offset for the 
pendulum model equal the initial center of mass offset for the CFD, but this approach did not 
necessarily improve the matching between the pendulum model and the CFD. If the number of 
pendulums in the model could vary and the model’s damping was non-linear, the approach of 
matching the initial center of mass offset may be more appropriate. 


A two-pendulum model is appropriate for all fill fractions. 


In some cases, a three-pendulum model that has the first pendulum dampen out quickly will 
give good results; in other cases, the behavior resembles a single pendulum model. In the case 
where a three-pendulum model seems more appropriate, the initial peak in forces is ignored 
and the two-pendulum model matches the subsequent peaks well. The 25% and 50% fill cases 
usually fell within this category. If a single pendulum is more appropriate, as was the case with 
85% fill, the frequencies were selected such that their constructive and destructive interference 
contributed to a better match of the CFD CM amplitudes. 


The forces and moments caused by the pendulum, as well as the equation describing the 
pendulum motion, as outlined in the Pendulum Model Derivation section. 


The force and moment equations, as well as the equation of motion for the pendulum, were 
derived from basic principles, using Dodge (2000) as a guide. This equation assumes that mass, 
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frequency, and damping are constant over time for a given propellant and fill fraction, although 
the CFD results indicate slight changes over time. The equation also assumes that damping is 
linear. 


The forces which cause the moment on the tank are located at the pendulum and not at the 
hinge point. 


Dodge (2000) suggests that the forces are applied at the pendulum location. 


CONCLUSIONS 


This paper developed a method for using two-pendulum models to predict slosh behavior on 
the Europa Clipper mission. The method resulted in a model that can be used to predict slosh 
behavior to a higher degree of accuracy than a single-pendulum model using less 
computational power than CFD. The parameters’ trends were supported by analytical and 
experimental literature. 


When implementing the pendulum models, it is crucial to recall its limitations. The model can 
accurately capture slosh behavior either before or after the first peak. High confidence is placed 
in the predicted pendulum frequencies, while moderate confidence is placed on mass, 
damping, and hinge locations. The model parameters may reflect inaccuracies in CFD, but 
experimental data is lacking to verify the CFD results. The pendulum model does not scale well 
for high fluid disturbance angles, but pendulum parameters based on a low disturbance angle 
will be conservative when scaled up to higher disturbance angles. The pendulum model 
assumes damping is constant over time, while it is in fact a function of time and distance 
traversed by moving fluid. The pendulum model is an approximation of true slosh behavior and 
these sources of error should be taken into account when designing a system to control slosh. 


The methods outlined in this paper could be further validated and refined if experimental data 
were available in low-gravity slosh experiments. 


CONTACT 
Wanyi Ng: wanyi.ng@nasa.gov 
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